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Abstract 

We investigate three Ising models on the simple cubic lattice by means of Monte Carlo methods and 
finite-size scaling. These models are the spin-|- Ising model with nearest-neighbor interactions, a spin-|- 
model with nearest-neighbor and third-neighbor interactions, and a spin-1 model with nearest-neighbor 
interactions. The results are in accurate agreement with the hypothesis of universality. Analysis of 
the finite-size scaling behavior reveals corrections beyond those caused by the leading irrelevant scaling 
field. We find that the correction-to-scaling amplitudes are strongly dependent on the introduction of 
further-neighbor interactions or a third spin state. In a spin-1 Ising model, these corrections appear to 
be very small. This is very helpful for the determination of the universal constants of the Ising model. 
The renormalization exponents of the Ising model are determined as yt = 1.587 (2), j/h = 2.4815 (15) 
and j/i = -0.82 (6). The universal ratio Q = (m^f/(m^) is equal to 0.6233 (4) for periodic systems with 
cubic symmetry. The critical point of the nearest-neig hbor spin-| model is IQ = 0.2216546 (10). 

Keywords: Phase transitions, Ismg model, Monte Carlo 
PACS: 64.60.Ak, 64.60.Cn, 64.60.Fr, 75.10.Hk, Ob.bO.+q 

1 Introduction 

According to insights such as the universality hypothesis, the nature of a phase transition does not depend on 
the microscopic details of a system but only on global properties such as dimensionality and symmetry of the 
order parameter. Thus, it is believed that most three-dimensional systems with short-range interactions and 
a scalar order parameter (such as density or unidirectional magnetization) belong to the Ising universality 
class. This implies that the critical exponents, as well as other universal quantities, are identical for all 
these models. This universality class comprises, in addition to anisotropic magnetic systems, also models for 
alloys, gas-liquid systems and liquid mixtures. 

In the case of two-dimensional Isinglike models, the evidence that universality holds is very strong. 
However, in three dimensions, where exact results are scarce and numerical techniques tend to be less accurate 
than in two dimensions, the situation is less satisfactory. Numerical uncertainties in the renormalization 
exponents amount to the order of several times 10~^. Since many years the most accurate results are those 
obtained by £-, coupling-constant and series expansions [1-9], whereas recently quite accurate estimates 
have been obtained by the coherent-anomaly method as well [10]. However, new possibilities to investigate 
Isinglike models are now arising in parallel with the availability of fast and relatively cheap computers. While 
many systems in the supposed Ising class may be simulated with the help of these, spin models offer a clear 
advantage, at least as far as a study of the universal properties is concerned. This is because of the ease 
and efficiency of the Monte Carlo method, in particular of cluster algorithms. Thus, results from Monte 
Carlo-based methods [11-20] tend to become increasingly accurate. 
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However, slight differences occur between recent results for the scaling dimensions. One possible expla- 
nation is that universality is not satisfied. In order to solve the issue whether these deviations are real, it is 
desirable to obtain more accurate Monte Carlo data for the supposed universal quantities. 

One problem that poses an obstacle to higher accuracies of these analyses is the presence of corrections to 
scaling. The dominant correction is attributed to an irrelevant renormalization exponent with an approximate 
value j/i ~ —0.83 [7]. This means that the corrections decay relatively slowly and thus jeopardize the accuracy 
of the analysis. For this reason, we explore what modifications of the simple cubic Ising model with nearest- 
neighbor interactions can infiuence the amplitude of these corrections to scaling. If we can, in this way, 
suppress the irrelevant field, we may expect a decrease of the ill effects due to the corrections to scaling. 
One can, for instance, choose a different lattice structure. Series expansions using the body-centered cubic 
lattice [7] indicate that corrections to scaling are relatively small. However, here we prefer to introduce 
continuously variable parameters to adjust the irrelevant scaling field. 

It is known [21] that the introduction of positive couplings with a range beyond the nearest neighbors 
in the simple cubic Ising model leads to a decrease of the correction-to-scaling amplitudes. We quote some 
preliminary results for the Hamiltonian 



where (nn) indicates a sum over nearest-neighbor pairs, (2n) over second-neighbor pairs (diagonals of the 
elementary faces), [3n] over third-neighbor pairs (body diagonals in the elementary cubes) and □ over four- 
spin products in all elementary faces of the cubic lattice. The associated couplings are denoted A'nn, K2n, 
A'sn and Kn, respectively. The spins Si can assume the values and —1. These results were obtained by 
Monte Carlo simulation on the Delft Ising System Processor [22,23]; they indicated that the introduction of 
positive A'2n, Ksn, or Kn reduces the correction-to-scaling amplitude. Third-neighbor couplings A'sn appear 
to be quite effective; for a ratio A'sn/A'nn ~ 0.4 the corrections become small. Much stronger second-neighbor 
couplings are required to obtain a similar effect [21]. 

Another approach is to introduce a third spin state Si = 0: the spin-1 Ising model. The weight of the 
Si = state can be varied by means of a term D s? in the Hamiltonian. Preliminary calculations showed 
that the corrections become small for D 0.7. In our actual simulations we have used _D = ln2, for reasons 
that will be explained in Sec. 2. 

Thus, we have selected the following three Ising models: the spin-i Ising model with A'sn/A'nn = (the 
nearest-neighbor model), with A'sn/A'nn = 0.4, and the spin-1 model with _D = In 2 and nearest-neighbor 
interactions. The algorithms used to simulate these models are described in Sec. 2. In addition to corrections 
to scaling, another obstacle to higher accuracies is the requirement of sufficiently accurate random numbers, 
in order to avoid biased results. In Sec. 3 we comment on the quality of our random-number generators 
and we mention consistency checks to which the algorithms were subjected. An analysis of the results for 
the dimensionless ratio Q = (m^)^/(m'*) is given in Sec. 4, followed by an analysis of the magnetic and 
temperature renormalization exponents in Sec. 5. The results for the three models satisfy universality: they 
are equal within the statistical inaccuracies. Assuming universality, Q as well as the critical points of the 
three models can be obtained with a better precision, as is demonstrated in Sec. 6. Finally, a discussion of 
these results in relation with the existing literature and with fundamental questions concerning universality 
is presented in Sec. 7. As an appendix, we present experimental results for critical exponents of a number 
of phase transitions that have been supposed to belong to the 3-d Ising universality class. 

2 Models and algorithms 

The present Monte Carlo analysis concerns three different Ising models. These can be represented in terms 
of a spin-1 Hamiltonian on the simple cubic lattice: 



Ti/kBT = - Ann ^ SiSj - K'zn ^ SiSj - Asn ^ SiSj - Kn ^ 
(nn) (2n) [3n] □ 



(1) 




(2) 
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Table 1: The ratio Ksn/ Knn and the value for D (see Eq. (2)) for the three models. 



model # 


A'sn/A'nn 


D 


description of model 


1 





— CO 


spin-i model with nn 


couplin 


gs 


2 


0.4 


— CO 


spin-i model with nn 


and 3n 


couplings 


3 





In 2 


spin-1 model with nn 


couplinj 





where nn and (ij) refer to nearest neighbors, and 3n and [kl] to third-nearest neighbors (along body diagonals 
of the elementary cubes). The spins can assume three discrete values Si = 0,±1. The three models are 
specified in Table 1. 

For D = — cxD the Si = states are excluded and thus models 1 and 2 can be simulated by the Swendsen- 
Wang (SW) [24], the largest-cluster (LC) [25] or the Wolff [26] method. In cluster algorithms, one has to 
'activate' a bond between two spins Si and Sj, coupled with strength Kij, with a probability p{Kij)bg^g^, 
where p(Kij) = [1 — exp(— 2A'jj)]. The presence of different sorts of bonds in model 2 thus leads to different 
bond probabilities but poses no further problems. If the bond is active, sites i and j belong to the same 
cluster. The simplest way to simulate this is to draw a random number for each bond and check whether it 
is smaller than p{Kij)bg^g^. Following this procedure, the speed of the algorithm decreases as the number of 
interacting neighbors increases. When the couplings are small, a more efficient procedure is possible. As a 
first step in the SW or LC cluster formation process one obtains, for each type of bond Kij, a list of bonds 
that should be activated if they connect equal spins. To this purpose, one introduces bond variables hij = 
or 1; the probability that hij = 1 is equal to p{Kij). The distribution P{k) = p(l — p)^~^, where we write p 
as an abbreviation for p(Kij), expresses the probability that (k — 1) subsequent bond variables equal zero, 
while the k-th bond variable is one. Thus one random number r can be transformed into an integer k: 

fc= l + [ln(r)/ln(l-p)] , (3) 

where the square brackets denote the integer part. After evaluation of k, the next (k — 1) entries in the list 
of bond variables are set to zero, and the k-th variable is set to one. By repetition of these steps a complete 
list of bond variables (for all bonds with strength Kij in the lattice) is obtained. Such lists are generated for 
each different type of bond. After completion of these lists, the cluster formation is trivial. This procedure 
was found to improve the speed of the simulation of model 2 considerably. One may still choose between the 
Swendsen-Wang or largest-cluster method. The latter method was observed to lead to shorter relaxation 
times and is therefore more efficient. The same principle was applied to Wolff-type simulations of model 2. 
Random numbers are, as above, transformed into integers k. During the cluster formation, (k — 1) bonds of 
the pertinent type are skipped and the spin connected to the k-th bond is added to the Wolff cluster if it has 
the right sign. This leads to a considerably faster Wolff algorithm, in particular because random-number 
generation is relatively time consuming (see Sec. 3). 

In the spin-1 case, transitions between zero and nonzero spin values require special attention. It is 
not immediately obvious how cluster algorithms could produce these transitions. We follow two different 
methods for the simulation of the spin-1 model. The first one uses a hybrid algorithm in which Metropolis 
sweeps alternate with cluster steps. The cluster algorithm acts on the nonzero spins only. Since we do not 
come close to the tricritical point where the ordered Ising phases meet the spin-zero phase, the regions of 
zero spins remain limited in size and we do not expect serious critical slowing down due to the equilibration 
between zero and nonzero spin values. 

The second method uses a mapping on a spin-^ model. We consider a Hamiltonian with two spins 
ti = ±1 and Ui = ±1 on site i (for all i) of the simple cubic lattice: 

-Hb/kBT = -Ml Y,(ii + + «i) - ^2 . (4) 

(ij) m 
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Using the transformation Si = (ti + Ui)/2 and Vi = (1 + ti)(l — Ui)lA, the partition function of this model is, 
up to a constant factor, 



where N denotes the number of spins in the system. This is, apart from the prefactor 2 , precisely the 
partition sum for Eq. (2) for A'nn = 4Mi, A'sn = and D = In 2 — 2M2. Equation (4) may thus serve 
for the application of cluster algorithms to the spin-1 Ising model. The special choice _D = In 2 leads to 
M2 = so that the spin-^ Hamiltonian simplifies. We have used three different methods to simulate the 
spin-1 model: the Metropolis-cluster (MLC) method, the full-cluster (EC) method and the Metropolis- Wolff 
(MW) method. The MLC method alternates one Metropolis sweep with one largest-cluster inversion, the 
MW method alternates one Metropolis sweep with 5 or 10 (this choice depends on the system size) Wolff 
steps. The EC method applies largest-cluster fiips to the spin-^ representation of the model: no Metropolis 
sweeps are included here. 

3 Random numbers and consistency tests 

Significant systematic errors may be introduced in Monte Carlo simulations by using inadequate random- 
number generators. It is well known that linear congruential methods based on the truncation of 32-bit 
integers are unsuitable for long simulations. Even their period of about 10® would be too restrictive. On the 
other hand, also random-number generators based on binary feedback shift registers may introduce serious 
errors (see, e.g., Refs. [22,27-30]). In most cases, the production rule selects two bits from the register 
and assigns their modulo-2 sum to the new bit. Thus the deviations from randomness are dominated by 
three-bit correlations. A number of algorithms of this type, using 127-bit shift registers with a period in the 
order of 10^®, have been rejected on the basis of long tests [23] using Metropolis simulations of the critical 
Ising model. Recent tests by Eerrenberg et al. [31] have shown that such deviations also occur when cluster 
algorithms are used together with random generators based on a generalized feedback shift register [32]. 

These findings apply to the Ising model on the square lattice, and use a comparison between simulations 
and exact results for finite systems. In the case of the three-dimensional Ising model, we have the practical 
difficulty that no exact results are available for general system sizes. One remaining possibility is a comparison 
between simulations with different realizations of the random-number generator and/or the spin-updating 
algorithm. Another possibility to check for systematic deviations is offered by numerical exact results for 
small systems. However, this test is not sufficient. It has been observed [33] that, in two dimensions, effects 
due to a random-number generator based on a generalized feedback shift register depend strongly on the 
system size and may be unobservable in small systems. 

Systematic effects in 3-d Ising model simulations are apparent from a comparison between finite-size 
results for the Binder cumulant [34] reported in Refs. [35] and [19], using the Swendsen-Wang and largest- 
cluster methods. The random generator used in Ref. [35] was based on a generalized feedback shift register 
with length 502 [36], and that used in Ref. [19] combined such a generator of length 9689 with a multiplicative 
rule, by means of bitwise modulo-2 addition. The observed discrepancies may be attributed to the shift- 
register-based generator with length 502 [36], and become even more prominent in simulations of the Wolff 
type [36]. 

It is clear that, for the long simulations implied by the present analysis of the 3-d Ising model, the 
random-number generators should be selected with great care. A systematic study of biases introduced 




(5) 



with Si = 0, ±1. Summation over the allowed values of yields a factor 2 if = 0. Thus 




(6) 
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by shift registers is necessary, in particular the dependence on the system size, shift-register length and 
the number of correlated bits. From simulations in two dimensions it appeared [37] that the deviations 
are scalable and become small for large system sizes and register lengths. Thus one may try to suppress 
systematic effects by using very long feedback shift registers [38]. But here we have chosen for a different 
method. This choice is based on the observation that the biases decrease when the number of bits in the 
production rule is increased [33,37]. The bitwise modulo-2 addition of two sequences generated by three-bit 
production rules usually leads to a sequence in which the dominant correlation is one between 9 bits. Thus, 
we expected that, using a random generator of this type with sufficiently long registers, the systematic 
effects would be well below the statistical accuracy, in three dimensions as well as in two. The largest 
part of the present simulations in three dimensions used the production rules Gi = ai-9218 © di-gesg and 
hi = &j_97 © &j_i27. These were combined by ri = Gi © hi, where Gi, hi and ri are 32-bit integers, and 
© stands for bitwise modulo-2 addition. Most of the simulations of the nearest-neighbor model reported 
in Ref. [19] were performed using a random generator which combines a multiplicative sequence with Gi. 
No systematic differences between both types of results were observed, nor were there obvious differences 
between simulations of the Swendsen-Wang, largest-cluster and Wolff types. Also in the case of models 
2 and 3 we checked for the presence of significant differences between the result of the different types of 
spin-updating algorithms (see Table 3) but none were found. This is consistent with the supposed high 
quality of the sequence ri . 

Therefore, we assume that the sequence ri is sufficiently uncorrelated, so that the simulation results 
may serve as a standard, to which data produced by means of other random generators can be compared. 
Thus, deviations in Wolff simulations of model 1, using 3-bit production rules, were determined and their 
scaling properties were analyzed [39]. The results are qualitatively the same as in two dimensions, and are 
completely consistent with the picture that the deviations decrease rapidly with increasing system size and 
shift-register length. No biases due to correlations of 5 or more bits were observed in Wolff simulations of 
the three-dimensional Ising model. 

As a further test, we have carried out exact numerical calculations of the dimensionless ratio QL(Knn), 
the susceptibility X; the energy-like quantity Snn and the specific-heatlike quantity Cnn- For more precise 
definitions we refer to the next sections. These calculations apply to small systems with periodic boundaries: 
model 1 with sizes 3^ and 4^, and model 2 with size 3^. 

The calculations for the 3^ systems involve a summation over 2^^ distinct states. In view of the efficiency 
of the calculation, we divided these states into subsets such that the states in each subset are related by 
symmetries: spin and spatial inversions, translations and rotations. The number of subsets is 55809, most 
of which contain 2592 states; this is the number of elements of the symmetry group of the 3^ Ising lattice. 
Once a list of 'independent' states, one per subset, has been made, the calculation becomes rather simple. 
However, the 4^ system is already too time-consuming unless treated carefully. The energy and the specific 
heat follow simply from the expansion coefficients given by Pearson [40]. For the magnetic quantities Q and 
X we have used a perturbation expansion similar to that described by Saleur and Derrida [41]. Details are 
given in Ref. [19]. The results are summarized in Table 2. The agreement between the Monte Carlo results 
and the exact numbers is quite satisfactory and does not suggest any problems with the random-number 
generator or other defects of the algorithm. 

4 Test of universality 

We have performed extensive simulations of models 1, 2 and 3, using the cluster methods described in Sec. 2. 
The total simulation time amounts to approximately two years on three workstations. We chose systems 
with size L x L x L and periodic boundaries. The lengths of the runs for the various models and methods 
are given in Table 3 for each system size. 

We sampled and analyzed the dimensionless ratio: 

04A-„„) = g| , (7) 
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Table 2: Comparison between Monte Carlo and exact results for small system sizes. These data were taken 
at couplings A'nn = 0.221653 and A'nn = 0.128006 for models 1 and 2 respectively. 



L 


model # 


quantity 


MC 


exact 


3 


1 


Q 


0.668409 (20) 


0.668427 


3 


1 


*> 

m 


0.422978 (19) 


0.422992 


3 


1 




1.434382 (51) 


1.434418 


3 


1 




0.785443 (53) 


0.785413 


4 


1 


Q 


0.659755 (24) 


0.659779 


4 


1 


*> 

m 


0.331228 (16) 


0.331204 


4 


1 


'5'nn 


1.293247 (40) 


1.293223 


4 


1 




0.977559 (75) 


0.977575 


3 


2 




0.642427 (20) 


0.642415 


3 


2 


*> 

m 


0.367413 (17) 


0.367390 


3 


2 


'5'nn 


1.134852 (51) 


1.134791 


3 


2 




0.287800 (19) 


0.287772 



Table 3: Length of Monte Carlo runs in millions of sampled configurations. SW stands for Swendsen-Wang, 
LC for largest cluster and W for Wolff. For SW and LC, each new configuration corresponds with one cluster 
decomposition of the lattice. In the Wolff case, 5 (5W) or 10 (lOW) Wolff clusters were fiipped before a new 
configuration was used for data taking. For the spin-1 model (model 3) LC, 5W and lOW are preceded by 
M in order to indicate a Metropolis sweep through the lattice. FC indicates the full-cluster algorithm for 
the spin-1 model; it fiips the largest cluster of a spin-^ version of the model. 



model 


1 


2 


3 


L 


SW 


LC 


5W 


low 


SW 


LC 


5W 


low 


FC 


MLC 


M5W 


MlOW 


3 


48 


52 


200 






100 


200 




100 


100 


300 




4 


160 


40 


200 






150 


100 




100 


100 


300 




5 


48 


52 


200 






150 


100 




100 


100 


400 




6 


48 


52 


200 






150 


100 




100 


100 


400 




7 


48 


52 


200 






150 


100 




100 


100 


400 




8 


48 


52 


200 




10 


140 


100 




100 


100 


400 




9 


48 


52 


200 




10 


140 


100 




100 


100 


300 




10 


48 


52 


200 




10 


140 


100 




100 


100 


300 




11 


48 


52 


200 






150 


100 




100 


100 


300 




12 


28 


72 


200 






150 


100 




100 


100 


200 




13 


28 


72 


200 






100 


100 








250 




14 


28 


72 


200 






100 


100 








200 




15 


20 


30 




200 




50 




100 






200 




16 


20 


30 




150 




50 




100 






200 




18 


12 


38 




150 




50 




50 






120 




20 


20 


10 




70 




50 




50 






120 




22 


18 


12 




70 




20 




80 






120 




24 


8 


12 




80 




20 




80 






150 




28 


10 


10 




100 




20 




80 






50 


50 


32 


2 


18 




180 




25 




75 








100 


40 




10 




90 




20 




80 








100 
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where L is the finite size of the model and m the magnetization density. We use the renormalization language 
in order to derive the expected finite-size scaling behavior ofQ^. By f(t, h, u, L~^) we denote the free-energy 
density as a function of the temperature and magnetic scaling fields, an irrelevant field and the finite-size 
field [42,43]. Here, we define the free energy as F = InZ, so without the normal factor — 1/A;bT. Its behavior 
under renormalization with a scale factor / is 

f{t, h, u, L-') = r^fiPH, P-h, P'u, l/L) + g{t, h) , (8) 

where j/t, J/h and m are the pertinent renormalization exponents, c? = 3 is the dimensionality and g is the 
analytic part of the transformation. By differentiating k times with respect to h, and choosing / = L and 
/} = 0, one obtains 

f'^^Xt, u, L-^) = L^y^-'^f'^^XLyH, L^'u, 1) + g^'^Ht) , (9) 

where the dependence on h is no longer needed and therefore suppressed. The expectation values of the 
second and fourth magnetization moments require differentiations of the free energy with respect to the 
physical magnetic field H: 

im') = L-'(U) (10) 

and 



4\ _ T-3d I 9 f \ , n r-^d f 9 f 



The Ising up-down symmetry implies that h is an odd function of H . Thus the correspondence between the 
derivatives with respect to h and H is 



and 



"^^-/"fl^V+V.^, ,13, 



where, as before, f^^^ stands for f / dh^ and all derivatives with respect to H are evaluated at iJ = 0. 
In the vicinity of the finite-size limit {t small and L finite), we may Taylor-expand the right-hand side of 
Eq. (9) in t and u. After the appropriate substitutions, the finite-size expansion of (5i(A'nn) follows as: 

QdK^n) = + ai(A'„„-A'c)i.^*+a2(A'„„-A'c)2i.2^*+a3(A'„„-A'c)^i.^^* + - • • + +62^^= + • • • , (14) 

where the Gi and hi are nonuniversal coefficients and t/2 = d—^^yh- The last term is due to the field dependence 
of the analytic part g in Eq. (9). The nonlinear dependence of h on H leads to even more rapidly decaying 
contributions (not shown). Terms of the same form, but with different exponents, may be due to other 
irrelevant fields. Because powers of the geometric factor dh/dH cancel in the first term, Q is a universal 
constant (related to the Binder cumulant [34]). 

The bulk of the numerical data were taken at couplings A'nn = 0.221653, 0.128006 and 0.393410 for models 
1, 2 and 3 respectively, close to the critical points. The results in terms of are shown in Table 4. A few 
points at somewhat different couplings were included in order to estimate the coefficients Gi in Eq. (14). The 
procedure of the analysis is as follows. We computed QL(Knn) for several values of _L, A'nn (near the critical 
points A'c) for the three models and fitted Eq. (14) to the data. The following parameters were used as 
input: j/t = 1.584 (4) (from e-expansion [5]; because the data were taken at couplings so close to the critical 
points, the results of the fits are practically independent of the precise value); j/i = —0.83 (5) (from series 
expansions [7]; the fit is rather sensitive to the precise value) and t/2 = —1.963 (3) (from renormalization 
arguments given above and the e-expansion result [5] for the magnetic exponent; the fit is insensitive to the 
precise value). The results are summarized in Table 5. It is stressed that the error margins quoted here 
include the uncertainty due to the possible variations in j/i, j/t and t/2 (Vh)- The fits for model 1 indicated that 
system sizes L < 7 should be discarded; they reveal finite-size effects not included in Eq. (14), exceeding 
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Table 4: Numerical results for the dimensionless ratio Ql = {rn^)\ I {rn^) l for the three Ising models defined 
in Sec. 2. These data were taken at couplings A'nn = 0.221653, 0.128006 and 0.393410 for models 1, 2 and 3, 
respectively. 



L 


model 1 


model 2 


model 3 


3 





66839 


(2) 





64244 


(3) 





61894 


(2) 


4 





65976 


(2) 





63164 


(3) 





62134 


(2) 


5 





65373 


(2) 





62642 







62242 


(2) 


6 





64919 


(2) 





62370 


f3) 





62273 


(2) 


7 





64579 


f3) 





62217 


f3) 





62277 


(2) 


8 





64318 


(3) 





62126 


(3) 





62288 


(2) 


9 





64107 


(3) 





62087 


(3) 





62280 


(3) 


10 





63943 


(3) 





62051 


(4) 





62280 


(3) 


11 





63803 


(3) 





62045 


(4) 





62275 


(3) 


12 





63688 


(3) 





62026 


(4) 





62272 


(3) 


13 





63591 


(4) 





62030 


(4) 





62267 


(5) 


14 





63514 


(4) 





62034 


(5) 





62262 


(5) 


15 





63441 


(3) 





62041 


(4) 





62266 


(5) 


16 





63376 


(4) 





62050 


(5) 





62247 


(5) 


18 





63270 


(4) 





62057 


(6) 





62248 


(7) 


20 





63187 


(6) 





62081 


(6) 





62230 


(8) 


22 





63117 


(6) 





62101 


(6) 





62211 


(8) 


24 





63052 


(6) 





62098 


(7) 





62201 


(7) 


28 





62958 


(6) 





62142 


(7) 





62177 


(9) 


32 





62879 


(5) 





62174 


(8) 





62129 


(8) 


40 





62761 


(8) 





62250 


(9) 





62050 


(9) 



Table 5: Results of a data analysis of the three models, including system sizes L > 7 for model 1 and L > 6 
for models 2 and 3. Besides the ratio Q, the critical couplings A'c and the nonuniversal coefficients ai, 02, 
hi and &2 are listed. 





model 1 


model 2 


model 3 


Q 


0.6232 (8) 


0.6229 (3) 


0.6231 (2) 




0.2216542 (8) 


0.1280034 (4) 


0.3934214 (8) 


ai 


0.862 (10) 


1.43 (4) 


0.659 (6) 


0-2 


0.54 (6) 


1.5 (2) 


0.352 (15) 


hi 


0.102 (10) 


-0.043 (4) 


0.001 (2) 


b2 


0.11 (3) 


0.351 (13) 


-0.018 (9) 
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the statistical error margins. The fits for models 2 and 3, which exhibit much smaller finite-size effects, 
include system sizes L > 6. The fit for model 2 clearly reveals a correction with exponent t/2 ~ —1.96. 
In fact, the large residuals in the absence of such a correction demonstrated its presence. As indicated 
above, this correction may arise from the analytic part of the transformation, although we cannot exclude 
contributions due to a second irrelevant exponent. Since there is no obvious reason why this term should 
be absent in general, we have included it in the fitting procedures for models 1 and 3 as well. Furthermore, 
we observe that the amplitude hi of the leading correction to scaling can be suppressed. This amplitude has 
become quite small in the spin-1 model (model 3) and has even changed sign in model 2. In model 1, the 
amplitude hi is relatively large and we have attempted to determine the irrelevant exponent by including it 
as a parameter in the fit. However, for an acceptable fit it was necessary to include the correction term h'^L^^ . 
Unfortunately, this frustrated the determination of m for model 1: if we fixed t/2 = —1.963 the exponent m 
shifted towards t/2 and if we included both t/i and t/2 as free parameters, they approached the same value. 

In order to take into account the finite-size effects revealed by the system sizes omitted in the previous 
fits, we have repeated our data analysis with an additional correction to scaling h^L^^ in Eq. (14), where 
t/3 = — 2t/h- This term, which is due to the nonlinear dependence of the magnetic scaling field on the 
physical magnetic field, enabled us to include all system sizes L > b for models 1, 2 and 3 in the analysis. 
The results, which are presented in Table 6, are consistent with those obtained previously. Again, the error 
margins quoted include the uncertainty due to the errors in t/i, t/t and t/h- These data satisfy universality 
within a margin of less than 10~^. To our knowledge, this is the most precise verification sofar for 3-d 
Isinglike models. 

Table 6: Results of a data analysis of the three models, where all system sizes L > b were used and a third 
correction term was included. Besides the ratio Q, the nonuniversal coefficients ai, 02, &i, &2 and 63 and the 
critical couplings A'c are listed. 





model 1 


model 2 


model 3 


Q 


0.6235 (7) 


0.6231 (4) 


0.6235 (3) 




0.2216547 (8) 


0.1280036 (5) 


0.3934224 (10) 


ai 


0.862 (9) 


1.43 (4) 


0.659 (6) 


0-2 


0.54 (6) 


1.5 (2) 


0.352 (15) 


hi 


0.098 (9) 


-0.045 (6) 


-0.004 (3) 




0.15 (4) 


0.37 (2) 


0.02 (2) 


&3 


-4.9 (8) 


-1.4 (8) 


-2.0 (7) 



5 Determination of the critical dimensions 

This section presents finite-size analyses of the energy, specific heat, spin-spin correlations over half the 
system size, susceptibility, the temperature derivative of the susceptibility, and the temperature derivative 
of the ratio Q^. Taking h = and / = _L in Eq. (8) leads to 

fit, u, L-i) = L-^fiiyH, Ly^u, 1) + git) . (15) 

Expansion in t and u yields 

fit, u, L-') = L-" (^/("■°) + f^''°^LyH + i/(2,o)|^2,,^2 ^ . . . ^ + f(^'^)Ly^+yHu + • • 

+g^°^+g^'h + lg(^h^ + -- - , (16) 

where /(*^'') stands for d^~^^ f / d^t d'u. The finite-size scaling behavior of the energy and that of the specific 
heat follow by differentiation. 
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5.1 The energy 

During the simulations, the nearest-neighb or sum Sun — ^!^(nn) ^i^j was sampled. For model 1, tliis sum is 
proportional to the energy; for models 2 and 3 its scaling behavior is similar. Its expectation value is equal 

(5„„) = ^ = ^^ + ^^. (17) 

The finite-size scaling behavior of this quantity thus follows by differentiating Eq. (16) and substitution in 
Eq. (17): 

(5nn) = Co + Ci(A'nn - A'c) H 

^ly^-d + a^(A'„„ - K,)Ly^ + a^{K^^ - K,fL^y* + ■■■ + + fes^^'"^* + • • •] , (18) 

where the Gi, hi and Ci are unknown coefficients. Analysis of the numerical results for (Snn) enables a 
determination of these coefficients and of j/f The dominant singular term in Eq. (18) is the one with 
amplitude oq. The (A'nn — A'c)-dependent term with amplitude ci is dominated by the term with coefficient 
ai and has therefore been omitted from the scaling formula. Since the bulk of the data were taken very close 
to the critical points, only linear and quadratic terms in (A'nn — A'c) were included. Without the correction 
term with coefficient &2, we had to exclude system sizes _L < 8 in the analysis of model 1, in order to obtain 
an acceptable residual. The resulting estimate for j/t is: 1.586 (6). Inclusion of the second irrelevant term 
enabled us to include all system sizes L > b. For consistency, we have included this term in the data analyses 
for models 2 and 3 as well. Table 7 summarizes the results obtained from fits according to Eq. (18), at the 
critical points listed in Table 6, for system sizes L > b. Since the singular behavior of (Snn) is rather weak, 
the results j/t ~ 1-59 for each of the three models are relatively inaccurate but consistent with the existing 
literature. The uncertainty due to the errors in A'c and j/i has been included in the error margins. 



Table 7: Results of a data analysis of the nearest-neighbor sum (Snn) for the three models. 





model 1 


model 2 


model 3 


yt 


1.599 (8) 


1.589 (9) 


1.591 (7) 


Co 


0.99051 (8) 


0.66298 (9) 


0.59451 (6) 


Go 


2.14 (6) 


2.20 (7) 


1.73 (4) 


hi 


0.14 (15) 


0.16 (17) 


0.04 (12) 


&2 


-2.0 (4) 


-0.6 (4) 


-0.9 (3) 



5.2 The specific heat 

The fluctuations in Snn are related to the speciflc-heatlike quantity 

nn 



A; 



5V 



™ (9A'2 

^ ni^ 

We consider / as a function of the scaling flelds t and u 

'df d'^t df d'^u 5V / dt 
'dt dK'L 5^ dK'L 5^2" VdK, 



d^f dt du d^f f du 



dtdu 9 A'nn 9 A'nn du"^ \dKn 



(19) 



(20) 



Taking the appropriate derivatives in Eq. (16) and collecting the leading analytic and singular terms leads 
to 

Cnn = P0+ Pi (A'nn " A'c) H 

+Ly^-d + _ K^'^LV. + ...]. (21) 
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The numerical results for Cnn of models 1-3 were subjected to a fit of this form with j/i = —0.83 and the A'c 
values in Table 6 as input parameters. The terms with amplitudes pi and si are dominated by that with 
amplitude qi and were omitted from the fit formula, as well as quadratic terms in (A'nn — Kc)- System sizes 
L < 6 display finite-size corrections not included in Eq. (21) and were discarded. The main results of these 
fits are shown in Table 8, where the error margins include the uncertainties in A'c and j/i. Also in the present 
case we find consistent, but inaccurate values of j/f This may be related to the fact that the leading power 
of L in Eq. (21) is close to zero, so that this term, which has the coefficient qo, interferes with the term with 
coefficient po- 

Table 8: Results of a data analysis of the specific-heatlike quantity Cnn obtained from the fiuctuations of the 
nearest-neighbor sum Snn for each of the three models. 





model 1 


model 2 


model 3 


yt 


1.60 (2) 


1.579 (15) 


1.59 (2) 


Po 


-0.8 (7) 


-0.6 (3) 


-3 (2) 


So 


1.5 (5) 


0.8 (3) 


3.5 (14) 


91 


2.2 (2) 


1.39 (12) 


3.7 (4) 


ri 


-0.4 (4) 


0.11 (14) 


0.0 (9) 


So 


-0.4 (3) 


-0.24 (13) 


-1.0 (7) 



5.3 The spin— spin correlation function 

In our simulations, we have sampled the spin-spin correlation function g(r), 

air) = (s(O)s(r-)) , (22) 

over half the system size (r = L/2), for even system sizes. This quantity can be derived from the free energy 
F by differentiating with respect to two physical magnetic fields Hq and Hr, which couple to the spins at 
positions and r, respectively. We consider the two fields as independent and find 

. d^F \ ^ f d^F dho dhr \ 

where h denotes the leading magnetic scaling field and the derivatives with respect to this field are evaluated 
at /lo = /ir = 0. Using Eq. (9) one obtains upon expansion in t and u the scaling behavior of the correlation 
function, 

g = L^y^-^" [ao + aiiK^n - K,)Ly^ + a2(A'nn - A'c)'!.'^* + • • • + + • • •] , (24) 

where the coefficients ai and hi are different from those in Eq. (18). 

We have fitted the terms shown in (24) to our data. The large residuals for all three models strongly 
suggested the presence of an additional correction to scaling h'^L^ . A problem for the determination of 
y' is the presence of the leading correction term Only in the spin-1 model (model 3), where the 

amplitude hi is small and the term thus may be omitted, a reasonable determination was possible, yielding 
y' = —2.1 (1). This could be a second temperature-like irrelevant exponent, although we haven't observed 
it in the analysis of the ratio Q or the energy-like quantity S-an- In Q, it may have been masked by the term 
h'^L^^ , but this is less likely for S-an , where the exponent of the second correction term is approximately equal 
to —2.4. On the other hand, the contribution h'^L^ could, in principle, be due to a second relevant magnetic 
exponent yh- Taking into account the dependence of F on an additional magnetic scaling field h yields 

52 A ( dhV f 52 A 52 A \ dh dh d'^F 



^^"■^ dhodK [dHj + [dhodK ^ dhodhj dHdH + dhodK ' ' ■ ^^^^ 
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This results in extra terms proportional to LV't^+y't^-'^'^ and i^2j/h-2d ^j^g scaling formula for g, corresponding 
to correction terms L^^~y^ and L^y^~^y^ in Eq. (24). Remarkably, the second magnetic exponent ijh = 
0.42 [14,18] has just the right value. However, its identification in terms of a redundant operator [14,18] 
would exclude its contribution to thermodynamic quantities. Table 9 shows the main results of an analysis 
for system sizes L > S, where we have included three correction terms, hiL^\ h'jL^ and h^L'^y , with the 
exponents m and y' fixed at —0.83 and —2.1, respectively. We have not included the term proportional to 
because hi is already quite small. The errors quoted in the table include the uncertainties in A'c, j/t, Hi 
and J/'. The estimates of j/h for each of the three models are consistent and in agreement with the existing 
literature. 



Table 9: Results of a data analysis of the spin-spin correlation function g for the three models. 





model 1 


model 2 


model 3 


J/h 


2.480 (2) 


2.482 (3) 


2.482 (3) 


ao 


0.77 (2) 


0.547 (15) 


0.453 (14) 


ai 


2.45 (4) 


3.07 (6) 


1.18 (2) 


0-2 


3.44 (14) 


7.2 (3) 


1.30 (3) 


hi 


-0.22 (10) 


0.08 (8) 


0.01 (7) 



5.4 The magnetic susceptibility 

The magnetic susceptibility x can be calculated from the average square magnetization, which is sampled 
in the Monte Carlo simulations, 

X = L^m") . (26) 
Using Eqs. (9), (10) and (12), we find for the finite-size scaling behavior: 

X = g^'\t) + L^y^-^'fi^HiyH, ly^u, i) , (27) 

which yields, upon expansion in t and u, 

X = co + ci(A'„„ - A'c) + • • • + L'y^-'' [ao + ai(A„„ - Ac)i.^* + a2(A'„„ - K,)L''y^ + + • • •] , (28) 

where the ai, hi and Ci are nonuniversal coefficients. In Table 10, we present the results of fits of the 
susceptibility for the models 1-3 at the critical points listed in Table 6. For model 1, system sizes L > 8 were 
included in the analysis and for models 2 and 3, which exhibit smaller corrections to scaling, all system sizes 
L > 6 were used. The coefficient ci in Eq. (28) was set to zero in all analyses, because the term containing 
it is much smaller than the (A'nn — A'c)-dependent term with amplitude ai. The errors include the margins 
due to the uncertainties in A'c, j/i and j/f The ratio between the coefficients ao for the three models is in 
excellent agreement with the ratio between the coefficients ao in Table 9. The same holds for the coefficients 
ai, a2 and hi. 

On the other hand, one might derive the scaling formula for x from that of the spin-spin correlation 
function g, because x is equal to the spatial integral of g, 

X = J g(ry-' dr . (29) 

Since the integral in Eq. (29) preserves the form of the corrections to scaling in g, we expect the same type of 
corrections in the correlation function and the susceptibility. Only the terms proportional to co, ci,. . . in (28), 
which arise from the analytical part of the free energy, are absent in Eq. (24). These contributions come 
from the small-r cutoff in Eq. (29). Thus, we have included in the scaling formula the additional corrections 
that we observed in the analysis of the correlation function. As the term proportional to interferes with 
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Table 10: Results of a data analysis of the susceptibility x for the three models. 





model 1 


model 2 


model 3 


J/h 


2.4812 (11) 


2.4817 (10) 


2.4826 (9) 


Co 


-0.6 (2) 


-0.20 (7) 


-0.50 (6) 


ao 


1.559 (16) 


1.126 (9) 


0.926 (7) 


ai 


4.88 (6) 


6.16 (8) 


2.36 (3) 


0-2 


6.9 (4) 


14.4 (7) 


2.62 (7) 


hi 


-0.37 (5) 


0.14 (3) 


-0.05 (2) 



the constant contribution cq, we have only included the correction h'zL"^^ . This allowed us to include system 
sizes _L > 5 for all three models. The results, which are presented in Table 11, are consistent with those 
obtained in the previous analysis. Now, the errors also include the margins due to the uncertainty in y' . Just 
as in the analysis of the correlation function, we find consistent results for j/h, which are in agreement with 
the literature. However, the values for j/h are more accurate than those obtained in the previous subsection 
and our resulting estimate for the magnetic renormalization exponent is j/h = 2.4815 (15). The error margin 
amounts to two standard deviations, in order to take into account any arbitrariness in the fit formula. 

Table 11: Results of a data analysis of the susceptibility x for the three models, where all system sizes _L > 5 
were employed and an additional correction to scaling was included in the scaling formula. 





model 1 


model 2 


model 3 


J/h 


2.4813 (11) 


2.4810 (14) 


2.4817 (13) 


Co 


-0.5 (2) 


0.0 (2) 


-0.33 (13) 


ao 


1.558 (15) 


1.134 (13) 


0.934 (10) 


ai 


4.87 (6) 


6.18 (8) 


2.37 (3) 


0-2 


6.9 (3) 


14.5 (7) 


2.64 (7) 


hi 


-0.37 (5) 


0.10 (6) 


0.01 (13) 


h2 


-5 (2) 


-2.9 (16) 


-2.6 (13) 



5.5 The temperature derivative of x 

In the simulations, we have also sampled the correlation between and Snn- This allows us to calculate 
the temperature derivative of the susceptibility, 

L'{{m'S^^)-{m'){S^^)) . (30) 



5A'„ 

The scaling behavior of this quantity can be derived directly from that of the susceptibility, Eq. (28), 
dx 



Cl 



+ . . . + i,2,,+,,-d I ^ ^ 2a2(A'„„ - K,)Ly' + 3a3(A'„„ - K,f 1""^' + • • • + liVJ' + • • • . (31) 



The term with amplitude hi comes from a term proportional to (A'nn — A'c)!^^*"*"^', included in the ellipsis 
in Eq. (28). Just as in the analysis of the spin-spin correlation function, the residuals for all three models 
indicated the presence of an additional correction to scaling h'^L^ , which indeed follows from the discussion 
in the previous subsection. Table 12 shows the results of an analysis at the critical points listed in Table 6, 
where this additional correction was included. All system sizes _L > 6 were used. The exponents m and y' 
were kept fixed at —0.83 and —2.1, respectively. The error margins include the uncertainties in A'c, j/i and 
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y' . The fit yields values for (2t/h + J/t) and the results for j/t have been obtained by fixing j/h at the best 
estimate from the previous subsection. This implies an additional error margin of 0.003 for j/f For models 
1 and 3, there is a reasonable agreement between the amplitudes ai and 02 as shown in Table 11 and those 
in Table 12. The differences are explained from the approximations in the scaling formulae. For model 2, no 
agreement is expected, because an additional term arises in the temperature derivative of the susceptibility 
due to the temperature dependence of A'sn (the ratio between A'nn and A'sn is fixed). 



Table 12: Results of a data analysis of the temperature derivative of the susceptibility dx/dKrm for the three 
models, where all system sizes L > 6 were employed and an additional correction to scaling was included in 
the scaling formula. 





model 1 


model 2 


model 3 


yt 


1.585 (3) 


1.584 (4) 


1.587 (4) 


Cl 


31 (11) 


1(9) 


12 (7) 


ai 


5.11 (9) 


3.73 (6) 


2.39 (5) 


0-2 


6.46 (10) 


8.04 (9) 


2.53 (3) 


0-3 


-2.4 (4) 


-4.4 (5) 


-0.53 (5) 


h 


-2.6 (3) 


0.1 (3) 


-0.1 (2) 


h 


-13 (3) 


-6 (2) 


-5.5 (16) 



5.6 The temperature derivative of Q 

Another quantity of interest correlates the magnetization distribution with the nearest-neighbor sum 5, 

dQ ^ ^ f Jm'Snn)-{m''){Snn) _ (m^g^^ ) - (m^) (g^^ ) 
dKnn \ {m?) {ni^) 

The determination of m and S-an during the simulations enables the sampling of this quantity with very 
little additional effort. Returning to Eq. (14) and noting that the ellipses include terms proportional to 
{K-an — A'c)!^^*"*"^' and to (A'nn — Kc)L^^ , we obtain the finite-size scaling behavior 

= [«o + «i(A'nn - K,)Ly^ + «2(Ann " Ac)'^'^* + • • • + t^i^' + + • • •] . (33) 

The numerical data for the three models were subjected to a fit on the basis of Eq. (33), where we have 
included system sizes L > 7 for model 1 and _L > 5 for models 2 and 3. In this case the leading power of L 
stands well apart from the less singular terms and the results for j/t (Table 13; uncertainties in A'c, j/i and 
j/h are included in all error margins) appear to be more accurate than those in the preceding subsections. 
The results suggest that the correction due to the leading irrelevant field is very small. Therefore we have 
repeated our analysis with v fixed to zero. We expect this to work especially well for models 2 and 3, 
where the irrelevant field is notably smaller than in the first model. Indeed, we have obtained accurate and 
consistent results for the models 2 and 3, as shown in Table 14. These results, together with those presented 
in Table 13, lead us to our final result, j/t = 1.587 (2). Just as in the final result for j/h in Sec. 5.4, we quote 
here an error margin of two standard deviations. 




6 Simultaneous fits for the three models 

Considering the results in the preceding sections, it is reasonable to assume now that universality is exactly 
satisfied for the three models under investigation. Thus we made a fit of the combined data for the ratio 
Q, allowing only single values of Q and j/i for the three models. The other parameters ai, 02, A'c, &i and 
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Table 13: Results of a data analysis of the derivative of the quantity Q with respect to the nearest-neighbor 
coupling A'nn • 





model 1 


model 2 


model 3 




1.589 (2) 


1.587 (2) 


1.5878 (14) 


Uo 


1.341 (14) 


1.351 (9) 


1.057 (6) 


Ml 


-0.21 (3) 


0.00 (5) 


-0.012 (9) 


M2 


-10.7 (6) 


-28.3 (11) 


-4.80 (12) 


V 


-0.01 (6) 


0.00 (3) 


0.02 (2) 


w 


-0.5 (2) 


0.46 (8) 


-0.13 (5) 



Table 14: Results of a data analysis of the derivative of the quantity Q with respect to the nearest-neighbor 
coupling A'nn for models 2 and 3, where the leading correction to scaling has been omitted. 





model 2 


model 3 


yt 


1.5868 (3) 


1.5867 (2) 


Uo 


1.3512 (11) 


1.0623 (7) 


Ml 


0.00 (5) 


-0.013 (9) 


M2 


-28.3 (11) 


-4.81 (12) 


W 


0.457 (15) 


-0.091 (9) 



&2 (see Eq. (14)) are nonuniversal and occur in triplicate. Now, system sizes L < 8 had to be discarded, 
except when an additional correction to scaling proportional to L~^^^ was added to the scaling formula. In 
the latter case, all system sizes _L > 5 could be included. Some of the results are summarized in Tables 15 
and 16, respectively, where the error margins include the uncertainty introduced by the error in j/t and t/h- 

Table 15: Results of a data analysis assuming universality of Ql = {rn^)\ I {rn^) l for the three investigated 
Ising models. System sizes L > 8 were included in the fit. The table lists nonuniversal parameters: the 
critical points and the amplitudes of the two correction terms. Furthermore, this analysis yielded the 
universal parameters Q = 0.6232 (2) and j/i = —0.78 (3). 





model 1 


model 2 


model 3 


Kc 


0.2216550 (6) 


0.1280037 (4) 


0.3934217 (8) 


hi 


0.086 (8) 


-0.040 (5) 


-0.001 (2) 


&2 


0.18 (3) 


0.34 (2) 


0.000 (14) 



Let us now compare the results of the various fits. In the first place, we see that the results in Tables 
15 and 16 are consistent, just as was the case for Tables 5 and 6 in Sec. 4. Also the values for the universal 
quantity Q, 0.6232 (2) and 0.6233 (2), respectively, agree. Secondly, the simultaneous fit with only the 
first two corrections to scaling (Table 15) yields results that are consistent with those presented in Table 5. 
Only the amplitude &2 and the critical coupling A'c for model 1 appear to be somewhat too low in Table 5, 
as we already had seen from the second fit in Sec. 4. Finally, when we compare the results in Tables 6 
and 16, i.e., including a third correction to scaling, as well as the corresponding Q values, we see a very good 
agreement. These comparisons, in addition to the fact that the term hsL^'^ allowed us to include all system 
sizes i > 5, lead us to the conclusion that the fits presented in Table 16 can be considered as the most 
accurate results. In addition to the nonuniversal constants given in the table and the universal amplitude 
ratio Q, this analysis yielded the (universal) irrelevant exponent m = —0.82 (3). This value is in very good 
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Table 16: Results of a data analysis assuming universality of Ql = {rn^)\ I {rn^) l for the three investigated 
Ising models. System sizes _L > 5 were included in the fit. The table lists nonuniversal parameters: the 
critical points and the amplitudes of the three correction terms. Furthermore, this analysis yielded the 
universal parameters Q = 0.6233 (2) and m = —0.82 (3). 





model 1 


model 2 


model 3 




0.2216546 (5) 


0.1280039 (4) 


0.3934220 (7) 


hi 


0.096 (7) 


-0.046 (5) 


-0.002 (2) 


&2 


0.15 (3) 


0.38 (2) 


0.007 (12) 


&3 


-4.9 (8) 


-1.7 (7) 


-1.7 (5) 



agreement with that obtained by Nickel and Rehr [7]. Although there is one more unknown (j/i), the results 
for Q and A'c obtained in this section are more accurate than those of the three separate fits. One of the 
reasons is that the fit for model 3 is insensitive to the value of m, so that, e.g., Q is determined accurately. 

7 Discussion and conclusion 

Let us summarize our final results for the renormalization exponents: j/t = 1.587 (2), j/h = 2.4815 (15), 
j/i = —0.82 (6). To allow for any residual dependencies on the choice of the fitting formulae, we list error 
margins of two standard deviations. In Table 17 we compare our results with some recent estimates obtained 
by various methods. Our result for the temperature exponent is lower than that of Baillie et al., obtained by 



Table 17: Some recent results for the renormalization exponents. 





Year 


yt 




Ui 


Present work 


1995 


1.587 (2) 


2.4815 (15) 


-0.82 (6) 


Kolesik and Suzuki [10] 


1995 


1.586 (4) 


2.482 (4) 




Guttman and Enting [9] 


1994 


1.580 (3) 






Landau [20] 


1994 


1.590 (2) 


2.482 (7) 




Baillie et al. [18] 


1992 


1.602 (5) 


2.4870 (15) 


-0.8 to -0.85 


Nickel [8] 


1991 


1.587 


2.4823 


-0.84 


Nickel and Rehr [7] 


1990 


1.587 (4) 


2.4821 (4) 


-0.83 (5) 


Le Guillou and Zinn-Justin [5] 


1987 


1.584 (4) 


2.4813 (13) 




Le Guillou and Zinn-Justin [2] 


1980 


1.587 (4) 


2.485 (2) 


-0.79 (3) 



the Monte Carlo renormalization method. This could be explained by a violation of hyperscaling. However, 
the accurate agreement between our result and that of coupling-constant expansion [2], e-expansion [5], 
series expansions [7, 8] and the coherent-anomaly method [10] makes this explanation less likely. We notice 
that our result for j/t is markedly higher than the recent series-expansion result of Guttman and Enting [9]. 
The result for the magnetic exponent is also in good agreement with most other estimates, although the 
result of Baillie et al. lies significantly higher than the majority of the results. Also the result of Le Guillou 
and Zinn-Justin obtained by coupling-constant expansion [2] seems somewhat too high. The results for the 
leading irrelevant exponent are not very accurate, but consistent. We notice that the fractions || and || are 
good approximations for j/t and j/h, respectively. For easy reference. Table 18 summarizes the exponents a, 
[3, 7, b, rj, V and A as calculated from our results for j/t, J/h and j/i, on the assumption that the hypotheses 
of scaling and hyperscaling are valid. 

Furthermore, we can calculate the Binder cumulant \J from our estimate for Q, using the relation \J = 
3 — l/Q, which yields \J = 1.3956 (10). Only few accurate results are available for this quantity (see, e.g.. 
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Table 18: The standard critical exponents as well as Wegner's correction-to-scaling exponent A as calculated 
from our best estimates for j/t, J/h and j/i. 



Exponent 


Expressed in RG exp. 


Value 


a 


2 - djy. 


0.110 (2) 


P 


{d- t/h)/t/t 


0.3267 (10) 


1 


(2t/h - d)/yt 


1.237 (2) 


6 


yh/{d-yh) 


4.786 (14) 


T] 


2 - 2t/h + 


0.037 (3) 


V 


1/j/t 


0.6301 (8) 


A 


-yjyt 


0.52 (4) 



Ref. [44] for a review) and one of the most accurate estimates up till now is U = 1.403 (7) [45]. Our result is 
in agreement with this and other estimates, but its accuracy is markedly higher. We have not sampled the 
characteristic length defined by Baker and Kawashima [46] in our simulations. Thus, we have no result 
for the renormalized coupling constant g* , which differs from the Binder cumulant by a factor {L/S^lY [46]. 

Table 19 presents a comparison of recent results for A'c of the spin-i nearest-neighbor Ising model. Again, 
it should be noted that the error margin of the result obtained in the present work amounts to two standard 
deviations. It can be seen that the amplitude ratio Q, which is used in this work, provides a good means of 
obtaining an accurate estimate for the critical coupling. We conclude that the conjecture of Rosengren [47] 
is not correct. The result of Ferrenberg and Landau deviates by 1.8 combined standard errors, but the 
newest estimate of Landau differs by only 1.2 standard deviations from the result presented here. The 
difference with Ref. [19] is 1.6 standard errors and is partly due to statistical errors (the data used in this 
work include those of Ref. [19] but are much more accurate and include _L = 40 data), and partly because a 
term with exponent t/2 was not included in the scaling formula for the ratio Q. Finally, we want to stress the 
importance of the spin-1 model. Since the corrections to scaling are small in this model, it is very suitable 
for the determination of universal quantities. 

Table 19: Summary of recent results for the critical point of the spin-^ Ising model with nearest-neighbor 
couplings. 



Reference 


Year 


Value 


Present work 


1995 


0.2216546 (10) 


Landau [20] 


1994 


0.2216576 (22) 


Blote and Kamieniarz [19] 


1993 


0.221648 (4) 


Baillie et al. [18] 


1992 


0.221652 (3) 


Livet [35] 


1991 


0.2216544 (10) 


Ferrenberg and Landau [17] 


1991 


0.2216595 (26) 


Ito and Suzuki [16] 


1991 


0.221657 (3) 


Blote et al. [13] 


1989 


0.221652 (5) 


Rosengren (conjecture) [47] 


1986 


0.2216586 (0) 
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A Experimental results 

In this appendix, we have collected a number of experimental results for the various critical exponents of 
phase transitions which have been compared to results for the 3-d Ising universality class. Only results 
published after 1980 are included. For binary mixtures, older results can be found in, e.g., Ref. [48]. The 
substances are grouped into four different subsets: unary systems, mixtures, magnetic systems and micellar 
systems (microemulsions). In general, the results agree very well with the theoretical values, but there are 
several remarkable discrepancies. Here, we only mention those measurements that differ by more than two 
standard deviations from our results. Results without error estimates are not taken into account. 

The value for 7 found in Ref. [69] lies much below the theoretical estimate. In Ref. [114] the results 
for both 7 and v do not appear to be Isinglike, as the authors have already noticed. The results for SFe 
presented in Ref. [62] are included, because the authors find mean-field values for the critical exponents, 
whereas the other results for the same substance appear to fall in the Ising universality class. For a, a range 
of values is found in Ref. [73], the lower end of which coincides with the theoretical value. In Ref. [104], two 
values of a are presented for MnF2, 2.4(7 above and 3.5(7 below the theoretical value, respectively. Several 
very accurate results for the critical exponent [3 are presented in Refs. [56,57,65], which all lie much above 
the result in Table 18. On the other hand, the accurate results in Refs. [94, 109] lie much lower than the 
theoretical prediction. The differences in Refs. [51,83,85] are less severe, but all results in these references 
lie between 2 and 3 combined standard errors above the theoretical value. Also the results for [3 presented in 
Ref. [95] deviate by many standard errors from the theoretical value, although the authors of this reference 
state that they consider the differences as not significant. The authors of Ref. [72] found a good fit of their 
data to a value of /3 = 0.3337 (5), but conclude that this value is probably too high. 

For v, very low values have been found in Refs. [84,93], whereas the value in Ref. [80] lies 2.6 standard 
errors above our result. In Ref. [118], the result rj = 0.016 (5) is given, which deviates by about four combined 
standard deviations from the best theoretical values. Nevertheless, the authors of Ref. [118] consider it to 
be in good agreement. 

In general, it is difficult to assess the source of the discrepancies noticed here, although there certainly 
are cases where crossover phenomena and corrections to scaling were not taken into account in the data 
analysis. 

Table 20: The various substances are grouped as follows: 'u' stands for unary systems, 'mx' for mix- 
tures, 'mg' for magnetic systems and 'mi' for micellar systems (microemulsions). a, [3, 7, b and v 
denote the standard critical exponents; A is Wegner's correction-to-scaling exponent; the critical expo- 
nents which carry a tilde refer to Fisher's 'renormalized' critical exponents [49]. Used abbreviations: 
AOT = di-2-ethylhexylsulfosuccinate; BHDC = benzyldimethyl-n-hexadecylammoniumchloride; CjEj = 
CH3-(CH2)i-r(0-CH2-CH2)j-OH; HFC-32 = difiuoromethane; HFC-125 = pentafiuoroethene; R114 = 
CC1F2-CC1F2; R13B1 = CBrFg; R12 = CCI2F2; R22 = CHCIF2. 



Substance 


Type 


Exponent 


Value 


Ref. 


CO2 


u 


a 


0.111 (1) 


[50] 








0.324 (2) 




Cs 


u 


a 


0.13 (3) 


[51] 






P 


0.355 (10) 




Ethane 


u 


P 


0.327 (2) 


[52] 






A 


0.46 (4) 




Fluoroform 


u 


P 


0.329 (1) 


[53] 


GeH4 


u 


P 


0.333 (8) 


[54] 


H2 


u 


P 


0.326 (3) 


[55] 
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oubstance 


Type 


Exponent 


Value 


Ret. 






7 


1.19 (5) 








A 


0.46 (2) 




HD 


u 


P 


0.352 (1) 


[56] 


HFC-32 


u 




0.345 (1) 


[57] 


HFC-125 


u 




0.341 (2) 


[57] 


Ne 


u 


P 


0.3575 (10) 


[56] 


Ne, N2 


u 


P 


0.327 (2) 


[58] 






A 


0.51 (3) 








A 


0.49 (5) 




R114 


u 


P 


0.312 


[59] 


R13B1 


u 


P 


0.340 


[60] 


R12 


u 


P 


0.337 


[61] 


R22 


u 


P 


0.348 


[61] 


Rb 


u 


a 


0.14 (3) 


[51] 






P 


0.36 (1) 




SFe 


u 


P 


0.48 (3) 


[62] 






P 


0.325 (5) 


[63] 






P 


0.338 


[64] 






P 


0.350 (4) 


[65] 






P 


0.355 


[66] 






1 


0.98 (5) 


[62] 






7 


1.24 (2) 


[65] 






b 


3.0 (2) 


[62] 


1,1, 1 ,2-tetrafluoroethane 


u 


P 


0.340 (1) 


[67] 


Xe 


u 


7 


1.246 (10) 


[68] 


Methanol + hexane 


mx 


7 


1.09 (3) 


[69] 






7 — a 


1.04(3) 




Methanol + cyclohexane 


mx 


P 


0.33 (2) 


[70] 






7 


1.26 (5) 


[71] 






V 


0.64 (2) 




Methanol + n-heptane 


mx 


P 


0.3337 (5) 


[72] 






a 


0.11-0.35 


[73] 


Methanol + isooctane 


mx 


P 


0.323 (9) 


[74] 


n-hexane + n-tetradecafluorohexane 


mx 


P 


0.35 (1) 


[75] 


Acetonitrile + cyclohexane 


mx 


P 


0.322 (4) 


[76] 


Butylcellosolve + H2O 


mx 


a 


0.077 (41) 


[77] 






P 


0.319 (14) 








7 


1.24 (1) 








V 


0.606 (18) 




Iso-butoxyethanol + H2O 


mx 


a 


0.105 (8) 


[78] 


Isobutyric acid + H2O 


mx 


P 


0.326 (3) 


[79] 






7 


1.19 (21) 


[80] 






V 


0.654 (9) 




Deuterated cyclohexane + cyclohexane + H2O 


mx 


P 


0.323 (3) 


[81] 






P 


0.326 (2) 








P 


0.322 (2) 




n,n-dimethylacetamide + octane 


mx 


P 


0.324 (5) 


[82] 






P 


0.329 (2) 




n,n-dimethylacetamide + decane 


mx 


P 


0.329 (4) 


[83] 






P 


0.333 (2) 
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oubstance 


Type 


Exponent 


Value 


Ket. 


Ethylammoniumnitrate + n-octanol 


mx 




0.610 (6) 


[84] 


Ethyleneglycolmonoisobutylether + H2O 


mx 


P 


0.332 (2) 


[85] 


2,6-lutidine + H2O 


mx 




0.336 (30) 


[86] 


Nitroethane + cyclohexane, benzonitrile + isooctane 


mx 




0.325 (5) 


[87] 


Nitrobenzene + isooctane 


mx 


a 


0.145 (35) 


[88] 


Nitrobenzene + decane + benzene 


mx 


P 


0.376 (8) 


[89] 


Perfluoroheptane + CCI4 


mx 


P 


0.324 (5) 


[90] 


Triethylamine + H2O 


mx 


A 


0.52 (3) 


[91] 


Triethylamine + H2O + D2O 


mx 


a 


0.110 (4) 


[92] 


Trimethylethylammoniumbromide + chloroform 


mx 


V 


0.621 (3) 


[93] 


Tetrachloromethane + tetradecafluoromethylcyclohexane 


mx 


P 


0.289 (6) 


[94] 


Tetra-n-pentylammoniumbromide + H2O 


mx 


P 


0.3370 (22) 


[95] 






P 


0.3190 (11) 








P 


0.3167 (16) 




CO2 + n-butane 


mx 


P 


0.359 


[96] 






V 


0.66 




CO2 + n-decane 


mx 


P 


0.368 


[97] 






V 


0.646 




Na + NH3 


mx 


p 


0.34 (1) 


[98] 






A 


0.46 (3) 




Na + ND3 


mx 


7 


1.228 (39) 


[99] 






7 


1.2400 (157) 








7 


1.223 (19) 








f] 


0.0300 (15) 








V 


0.0317 (13) 


[100] 






V 


0.0302 (15) 








v 


0.6279 (80) 


[99] 


C0F2 


mg 


a 


0.109 (6) 


[101] 


DyAlG 


mg 


P 


0.33 (1) 


[102] 


FeCl2 


mg 


a 


0.15(4) 


[103] 


FeF2 


mg 


a 


0.111 (7) 


[104] 






a 


0.115 (4) 








a 


0.11 (3) 


[105] 






P 


0.325 (2) 


[106] 






7 


1.25 (1) 


[107] 


MnBr2 


mg 


a 


0.118 (7) 


[108] 


MnCl2 


mg 


P 


0.297 (3) 


[109] 


MnF2 


mg 


a 


0.123 (5) 


[104] 






a 


0.091 (5) 




NdRu2Si2 


mg 


a 


0.11 (3) 


[110] 


U3P4 


mg 


P 


0.315 (15) 


[111] 






P 


0.313 (15) 








7 


1.25 (2) 




AOT + n-decane + H2O 


mi 


7 


1.26 (10) 


[112] 






7 


1.22 (5) 


[113] 






7 


1.61 (9) 


[114] 






v 


0.61 (6) 


[112] 






V 


0.75 (5) 


[113] 






V 


0.68 (8) 








V 


0.72 (4) 


[114] 



20 



PREPRINT 



fsuDstance 


Type 


Exponent 


Value 


itei. 


Benzene + BHDC + H2O 


mi 


/? 


0.34 (8) 


[115] 






7 


1.18 (3) 


[116] 






v 


0.60 (2) 


[116] 


2-butoxyethanol + D2O 


mi 


1 


1.216 (13) 


[117] 






V 


0.039 (4) 








v 


0.623 (13) 




CeEs + H2O 


mi 


a 


0.11 (4) 


[118] 






P 


0.327 (4) 








1 


1.24 (1) 








1 


1.241 (16) 


[119] 






V 


0.016 (5) 


[118] 






v 


0.627 (6) 








V 


0.632 (11) 


[119] 


C8E4 + H2O 


mi 


1 


1.237 (7) 


[119] 






1 


1.243 (7) 








V 


0.630 (12) 








V 


0.630 (18) 




C10E4 + H2O 


mi 


1 


1.25 (2) 


[120] 






V 


0.63 (1) 




C12E5 + H2O 


mi 


1 


1.17 (11) 


[121] 






V 


0.65 (4) 




C12E6 + D2O 


mi 


1 


1.2(1) 


[122] 






V 


0.60 (5) 




C12E6 + H2O 


mi 


1 


1.2(1) 


[123] 






V 


0.60 (3) 




C12E8 + D2O 


mi 


1 


1.21 (2) 


[124] 






V 


0.62 (2) 




C12E8 + H2O 


mi 


1 


1.20 (4) 


[124] 






V 


0.63 (1) 




Sodiumdodecylsulfate + butanol + NaCl 


mi 


V 


0.62 (3) 


[125] 






V 


0.64 (4) 








V 


0.63 (5) 




Cationic surfactant in aqueous salt solution 


mi 




0.375 (10) 


[126] 






1 


1.39 (4) 








V 


0.70 (3) 
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